Anomalous heat conduction in one dimensional momentum-conserving systems. 



(N 
O 

o 

(N 



Onuttom Narayan 1 ' 2 and Sriram Ramaswamy 1 
1 Centre for Condensed Matter Theory, Department of Physics, 
Indian Institute of Science, Bangalore 560012, INDIA 
2 Department of Physics, University of California, Santa Cruz, CA 95064-* 
(Dated: February 1, 2008) 

We show that for one dimensional systems with momentum conservation, the thermal conductivity 
k generically diverges with system size as L 1 ^ 3 . 
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When a very small temperature difference is applied 
across a system, it is expected that in steady state the 
heat current j will obey Fourier's law of conduction 



3 = -«VT 



(1) 



where T is the local temperature and k is the heat con- 
ductivity of the material. Although k is in general tem- 
perature dependent, if the applied temperature difference 
is small it should be constant across the system. Thus if 
Xi and T 2 are the temperatures at which the two ends of 
a system of length L are kept (with T\ > T 2 ), the steady 
state current should be j = n{Ti — T 2 )/L. 

On the other hand, for many one-dimensional mod- 
els it can be shown analytically f5L @, || or numeri- 
cally H [| [| |, § that j cx L a - X with a > 0, even in the 
linear response regime where j oc Ti — T 2 . The value of a 
differs from model to model. In the framework of Eq. ([!]) , 
this would imply an L-dcpcndent conductivity that di- 
verges in the infinite system limit. (In some oscillator 
models, a < 0, implying an anomalous but not divergent 
conductivity.) 

Recently, it has been argued || that such anomalous 
heat conduction occurs only in systems with momentum 
conservation, and is a consequence thereof. This was 
done by showing that if j(x, t) is the energy current den- 
sity, the autocorrelation function of the total energy cur- 
rent J(t) = J dx j(x,t) 



C(h-t 2 ) = {J(t 1 )J(t 2 )) 



(2) 



has the property C(t — > 00) ^ 0. Although the proof 
of this result in Ref. || was specific to one dimensional 
systems, Galilean invariance allows one to construct a 
general proof for any dimension d This is because 
the energy current J(t) has an advective contribution 
(E + pL d )v, where E is the energy, p the pressure, L d 
the volume and v the center of mass velocity of the 
system. For an energy and momentum conserving sys- 
tem, this advective contribution to J(t) is time inde- 
pendent. C(oo) is found by calculating ([(E + pL d )v] 2 ) 
within the canonical ensemble, and is non-zero. The 
limit lim T ^oo liiri£_, 00 (l/ (T 2 L)) J Q dtC(t) is then diver- 
gent. The Kubo formula was invoked || to equate 
this to k. 

Although (as discussed in the next paragraph) the ar- 
gument in Ref. @] is incorrect and a nonvanishing C(oo) 



has no consequence for heat conduction, the conclusion 
that momentum conservation in low dimensional sys- 
tems (generically) implies anomalous conduction is valid. 
This is because, in addition to the limit C(co) being 
finite, there is also a slowly decaying tail in Cb(t) = 
C(t) — C(oo). Unlike C(oc) 7^ 0, which is valid in all 
dimensions, the tail in Co(t) decays sufficiently slowly to 
cause a singularity in k only for d < 2. In this paper, 
using the transport equations for a normal fluid with 
thermal noise added, we show that k(L) oc L a , with 
a = (2 - d)/(2 + d) > for d < 2. (There is a loga- 
rithmic singularity for d = 2.) For the physical case of 
d = 1, we obtain a = 1/3. These transport equations 
should be valid for all systems which rapidly reach lo- 
cal thermal equilibrium, and for which the only slowly 
evolving quantities are the mass energy and momentum 
densities. The results are therefore generic, although the 
assumption of local thermal equilibrium breaks down for 
some models, such as hard sphere equal mass particles 
in one dimension, or a chain of harmonic oscillators, for 
which a is different §. 

We first recall the discussion of Bonetto et al Jll[] 
about the argument in Ref. p|. The conductivity in 
momentum conserving systems is in fact not obtained 
from the autocorrelation function of J(t) , but of J (t) = 
J(t) — (E + pL d )v. |1| This result can be proved rigor- 
ously within linear response theory |l3| . If one makes a 
linearized hydrodynamic approximation, which is equiv- 
alent to linear response theory with no (thermal) noise 
terms in the transport equations, the derivation of this 
result is simpler and well known jlj, p^| . When J(t) 
is replaced by Jo(t), the resultant truncated correlation 
function Co(t) has no infinite time tail, and 



k = lim lim — = — 

T^OO L-^OO T L 



dtC (t) 



(3) 



is cured of its divergence. 

Despite the fact that Co (00) = 0, we show in this 
paper that, with thermal noise in the transport equa- 
tions, Co(t) decays sufficiently slowly with t to make 
the integral J dtCo(t) — and therefore the conductiv- 
ity k — divergent for d < 2. We first show this qualita- 
tively. The advective contribution to J(t) is really equal 
to J d d xh{x)v{x), where h(x) is the local enthalpy den- 
sity and v(x) the local velocity. (Vector indices have been 



2 



suppressed.) This can be expressed as J d d kh(k)v(—k) + 
(E+pL d )v, where the first integral is restricted to k ^ 0, 
and the second part comes from the non-zero spatial 
average of h{x). In going from J(t) to Jo(t), only the 
second part of this was removed. Since the time decay 
of all the hydrodynamic modes is diffusive p^| , express- 
ing h(k,t) and v(k,t) in terms of hydrodynamic modes 
and approximating (h(k, t)v(— k, t)h(k', 0)v(— k', 0)} as 
(h(k,t)h(k',Q)){v(-k,t)v(-k',0)) yields the advective 
contribution to C (t) to be ~ L d j d d k exp[-0(k 2 )t], 
which is ~ L d /t d / 2 . From Eq.(|), k diverges for d < 

Although the conductivity n does indeed diverge for 
d < 2 as this rough calculation indicates, the tail of Co(t) 
decays as £-2d/(<i+2) instead of as t~ d / 2 . This is because 
thermal noise in the transport equations gives rise to sin- 
gular corrections to the parameters in the equations, so 
that the hydrodynamic modes decay superdiffusively. As 
is standard in renormalization group (RG) analyses of 
such phenomena, we solve the linearized transport equa- 
tions and check whether the non-linear corrections are 
relevant for long wavelength low frequency phenomena. 
Below d — 2, which is thus the upper critical dimension, 
the nonlinearities are found to be relevant, and their ef- 
fect is calculated. 

We assume that the system whose thermal conductiv- 
ity is of interest is one which reaches local thermal equi- 
librium, with the only dynamical variables that evolve 
slowly with time being the mass, energy and momen- 
tum densities. With these assumptions, the appropriate 
transport equations for the system are those for a normal 
fluid Jl7| , with thermal fluctuations included in the form 
of noise sources [[l9| |o| : 

d t p + w-{pv) = 

d t (pv) + V • (pvv) = -Vp + (C + ?7/3)VV • v 
+ rjV 2 v + Cv 

d t e + V- [{e + p)v] = V-k VT + O((Vi;) 2 ) + C £ (4) 

where p is the local density of the fluid. The lo- 
cal temperature T and pressure p are implicit func- 
tions of p, v and e. The thermal noise terms £„,e satisfy 
((v(xi,h)(v(x 2l t 2 )) oc -k B T8(ti - t 2 )d 2 5(xi - x 2 ), and 
similarly for £ c , with the proportionality constants fixed 
by the requirement that the variance of the fluctuations 
bp = p — po, 5e = e — €q and v at any instant are those of 
a system in equilibrium at temperature T. (V • (pvv) is a 
vector whose i'th component is dj(pViVj). The first equa- 
tion in Eqs.(^j) is an exact identity, and has no thermal 
noise correction. 

It is standard and straightforward to solve these equa- 
tions with the linear approximation, valid if Sp(x, t), 
v(x,t) and Se(x,t) are small. One obtains [ fbl] , [l7f two 
propagating sound modes and one non-propagating heat 
mode. All three modes decay diffusively. The relevance 
or irrelevance of the nonlinear terms in the transport 
equations is then determined by calculating corrections 
to correlation and response functions to one loop in a 



diagrammatic perturbation expansion. For instance, the 
response of v(x) to a perturbation in the second mem- 
ber of Eqs.(^), g vv (x\ — x 2 ;ti — t 2 ), receives a one loop 
correction from V • (pvv), resulting in corrections to the 
viscosities 77, £ of the form 

5r],8£~ J d d xdtc vv (x;t)g vv (x;t) (5) 

where c vv is the autocorrelation function of v(x). (The 
component indices in c Vi ^ . and g Vi }V . have been sup- 
pressed for compactness.) Expanding the correlation 
function c and the response function g in terms of the 
three hydrodynamic modes, and performing the x inte- 
gral first, the integrand is negligible outside a region of 
volume 0(t) d / 2 . (The propagating parts of the modes 
shift the peak of the integrand away from x = if the 
contribution of the same hydrodynamic mode is consid- 
ered for c and g. This is inconsequential if the system is 
large.) Since the correlation and response functions of 
all the hydrodynamic modes have a \t\~ d / 2 prefactor, the 
integral over x yields Srj, <5£ ~ J dt9(t)t~ d / 2 , where the 9- 
function comes from causality in the response function. 
The t integral diverges for d < 2. Similar calculations can 
be carried out for other one loop corrections. Since all 
the autocorrelation functions have the same \t\- d / 2 pref- 
actor, the RG scaling dimension of all the three density 
fields is -d/2 : \t\- d / 2 - (\x\- d / 2 ) 2 . 

For d < 2, the nonlinear corrections are therefore 
relevant when expanding around the linearized equa- 
tions. This can also be seen by scaling all the vari- 
ables in the transport equations as x — Xx' , t — X z t', 
C„,e = A-( d+2+z )/ 2 C; e , and (Sp,v,Se) = \- d / 2 (Sp,v,Se). 
The time derivative, dissipative and thermal noise terms 
in Eqs.(||) scale identically if the dynamic exponent z is 
set to 2. (The terms with one spatial derivative in the 
linearized equations grow, since they control the prop- 
agation of the sound modes, whose speed is obviously 
altered if t is scaled as x 2 . However, as we have seen in 
the previous paragraph, this does not affect the scaling 
of the loop corrections.) The nonlinear terms can be seen 
to be relevant if d < 2. 

For d < 2, a renormalization group analysis has to be 
carried out, integrating out loop corrections from short 
wavelength fluctuations along with rescaling the vari- 
ables. The nonlinearities grow under the RG flows un- 
til they reach a non-trivial fixed point. The new scal- 
ing dimensions of the fields and the dynamic exponent 
z can then be evaluated at this fixed point. It is, how- 
ever, not necessary to carry out such a calculation to ob- 
tain the exponents: they can be determined completely 
from symmetry considerations. The RG flows preserve 
the property that equal time fluctuations in 8p, v and Se 
must be those of a system in equilibrium at tempera- 
ture T. Since the fluctuations in these densities must be 
Gaussian at sufficiently long wavelength, we require that 
J d d x[v 2 , (Sp) 2 , (Se) 2 } should be invariant under rescaling. 
Thus the scaling dimensions of all three fields are equal 
to —d/2 even for d < 2. Further, Galilean invariance re- 



3 



lates the loop corrections to dt4> and to the corresponding 
advective term V • (v<p) for any (conserved) field <f> |2lf| 
Since both terms are invariant at the fixed point, this 
yields the condition that Vv scales as dt, i.e. v scales as 
x/t. Combining this condition with the previous one, we 
obtain z = 1 + d/2. 

The energy current density is obtained by requiring 
that the third equation in Eq. (^J) should be equivalent to 
d t e + V • j = 0. This yields j(x, t) = (e + p)v - k VT + 
0(Vv 2 , wV • v). Since J (t) = J(t) - {E+pL d )v, we ob- 
tain correspondingly jo(x, t) = (Se + p — po)v — k VT + 
0(Vv 2 ,vV ■ v). Under the RG rescaling, the three terms 
in this scale with dimension — d, — 1 — d/2 and —1 — d 
respectively, and the first term is most important for 
d < 2. If one expresses the transport equations Eqs.(§) 
through a generating functional p2] , and adds an extra 
a(Se + 5p)v term in the argument of the exponent, the 
v — ► — v,x — ► —x symmetry of Eqs.(^) ensures that this 
term is not renormalized to 0(a). Thus the scaling of 
Co(t)/L d can be obtained from the bare scaling dimen- 
sion of jo(x, t) : 

C (t)/L d = J d d x{j o (x,t)j o (0,Q)) ~ Itl"- 1 (6) 

with 

a=l-d/z=(2-d)/(2 + d). (7) 

Eq.(||) has to be integrated over t to obtain k from 
Eq.(0). For a system of linear dimension L, the tail of 
Co (i) obtained in Eq. (^|) is valid only when a disturbance 
at x = 0, t = in the propagating modes has not been 
carried outside the system. This is because the tail in 
(h(x,i)v(x,i)h(0,0)v(Q,0)) comes from long wavelength 
fluctuations, for which the contribution to v from the 
heat diffusion mode is zero, so that v must be expressed 
as a linear combination of the two propagating modes. A 
fluctuation in v is carried to the boundaries of the sys- 
tem in a time t ~ O(L). The precise behavior thereafter 
depends on the coupling to the heat baths at the bound- 
aries, but in any case the fluctuation is partially or fully 
lost to the baths. The tail of Co(t) is cut off in a few 
round trip times, i.e in a time t ~ O(L). Substituting 
Eq.(§) in Eq.(§) and using this cutoff, we obtain 

k(L) = L a . (8) 

Thus the conductivity measured in a system of size L 
diverges with L, or equivalently, the heat current flowing 
across a system with a fixed small temperature difference 
decays as ~ with L. For the physically relevant case 

of d = 1, the heat current must decay as ~ L -2 / 3 . For 
d = 2, the conductivity n has a logarithmic singularity 
as a function of L. 

We now compare with earlier analytical and numeri- 
cal results. For a chain of coupled harmonic oscillators, 
there has been a large amount of analytical work show- 
ing that k diverges with system size [[l] . The form of the 
divergence is different for different models, and in fact 



varies over a wide range depending on the heat baths M . 
However, all these are systems where local thermal equi- 
librium is not established, so the results of this paper do 
not apply. With more complicated models, there have 
been various numerical simulations that have shown a 
divergent k, with a ranging from 0.17 to 0.5 [|[ ||, [|. 
However, the simulations in these papers were performed 
for fairly small system sizes, up to ~ 1000. The scaling of 
k as a function of L is not very good, indicating the need 
for larger system sizes. Very recently, there have been 
simulations on a one dimensional hard sphere gas with 
alternating masses p3| for much larger system sizes: up 
to 16383 § and 30000 @]. In the former, the scaling is 
not very good, and the dependence of k on L varies con- 
siderably with the mass ratio between neighboring parti- 
cles. However, the authors estimate a to be 0.31 to 0.35. 
In the latter paper, the scaling is very good, and a is 
0.255, which disagrees with our paper. It is surprising 
that numerics on the same model, with roughly the same 
range of sizes, yields such different results. The current 
j{x,t) is chosen differently in both papers, and it is not 
clear how it is defined in the latter paper, since it in- 
volves derivatives of the (singular) hard sphere potential. 
Further work is needed to clarify the situation. 

In Ref. [|, there is a mode-coupling calculation that 
relies on Ref. |24|j , indicating that a should be 2/5. The 
argument in |24| is internally inconsistent: the scaling 
u) ~ D(u})k 2 is used to correctly find the renormalized 
diffusion coefficient D{uj) ~ w -1 / 3 and thereby the long- 
time tail of Co(t) as ~ i~ 2 / 3 , but then w ~ k is used to 
incorrectly convert D{u>) to D{k) ~ fc -1 / 3 . This expres- 
sion is then used by Ref. J| to obtain C (t) - £~ 3/5 . As 
we have seen, although t ~ x is the scaling conversion to 
be used in Eq. (||) , this is not appropriate for the loop inte- 
grals or the dynamic exponent. Secondly, although their 
system is nominally a 1-dimensional crystal, it should be- 
have like a fluid at large length scales since fluctuations 
will wipe put long-range order. Nonetheless, if we treat 
it as a crystal, we will encounter the "Poisson-bracket" 
nonlinearity S/uSH/Su in the equation for the momen- 
tum density, where H is the elastic Hamiltonian for the 
displacement field u. This term is nonlinear even if we 
retain only harmonic terms in H. Simple power-counting 
shows that this yields nonlinear corrections of precisely 
the same type as those from the advective term. In fact, 
the analysis of Ref. || is equivalent to this. 

For a system in which momentum is not conserved, 
there is no advective term in the energy or mass current 
both of which depend on gradients of p and e. Even the 
nonlinear terms in the equations for dtp and dte thus 
have at least two spatial derivatives, and are irrelevant 
compared to the linear terms. The conductivity is given 
by Eq.(|) with C{t) instead of C (t) || || This neither 
has a non-zero C(oo) limit, nor a slowly decaying tail for 
large t. Indeed, numerical studies pf| confirm that for 
such systems, the conductivity is finite. 

If a system is integrable, even if it does not con- 
serve momentum, the conductivity can be singular, since 
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Refs.|l3|, assume local thermal equilibrium. Also, re- 
cent numerical studies on a momentum conserving chain 
of coupled rotators |2rj| show no anomalous conductiv- 
ity at high temperatures. The interparticle potential is 
V(gi+i — qt) = 1 — cos(<7i + i — qi), and the particles at the 
end are kept in contact with different heat baths (through 
Langevin noise). Physically, the q's can be interpreted 
as angles. Although the model is momentum conserv- 
ing, since the heat baths are applied to specific particles 
rather than at the ends, there is no advective term in the 
current: a long wavelength momentum fluctuation causes 
the "angles" to "spin" round and round instead of car- 
rying energy from one side to another. At low temper- 
atures, when fluctuations are small and the model goes 
over to a harmonic oscillator, a long wavelength momen- 
tum fluctuation is transmitted to neighboring particles 



and thence to the ends of the chain, effectively leading 
to advective transport. The nature of the transition be- 
tween the high and low temperature phases is interesting 
and requires further work. Also, if a system is integrable, 
even if it does not conserve momentum, the conductivity 
can be singular, since Rcfs.Jl^, [l^] assume local thermal 
equilibrium. 

In this paper, we have shown that for one-dimensional 
momentum conserving systems, the heat current when a 
small temperature difference ST is applied across a sys- 
tem of length L is generically of the form j oc (ST) / L 2 ^ 3 . 
This is consistent with earlier numerical studies, but fur- 
ther work is needed to improve the numerical picture. 

We thank Abhishek Dhar and Sriram Shastry for very 
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